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Abstract 

We study the elasticity of perfect ^He at zero-temperature using the diffusion Monte Carlo 
method and a realistic semi-empirical pairwise potential to describe the He-He interactions. In 
particular, we calculate the value of the elastic constants of hep helium {Cij} as a function of 
pressure up to ~ 110 bar. It is found that the pressure dependence of all five non-zero {Cij} is 
linear and we provide an accurate parametrization of each of them. Our elastic constants results are 
compared to previous variational calculations and low-temperature measurements and in general 
we find notably good agreement among them. Furthermore, we report T = results for the 
Griineisen parameters, sound velocities and Debye temperature of hep ^He. This work represents 
the first of a series of computational studies aimed at thoroughly characterizing the response of 
solid helium to external stress-strain. 

PACS numbers: 67.80.-s,02.70.Ss,67.40.-w 
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I. INTRODUCTION 



The behavior of sohd ^He, and of quantum crystals in general (e.g., H2 and Ne), is 
exceptionally so rich that despite having been investigated for more than about eight decades 
is to this day not yet completely understood. One example of helium's intriguing nature 
is its elasticity. Experimental studies on the elastic properties of hep ^He were already 
conducted by Wanner, Crepeau and Greywall in the early seventies.-"- Those original works 
consisted of series of sound-velocity measurements performed at thermodynamic conditions 
relatively close to the stability domain of the liquid, namely T ~ 1 K and 25 < P < 50 bar. 
With the advance of time and technology cryogenic and crystal growth techniques have been 
improved so notably that at present is possible to analyze practically defect-free ^He samples 
at just few mK in the laboratory. Recently, Beamish and collaborators have developed a 
new experimental technique that has allowed them to measure directly the shear modulus 
fi of hep ^He under extremely low strains and frequencies.™!^ The temperature dependence 
of fi within the temperature interval 0.5 < T < 0.01 K has been determined and a striking 
resemblance with non-classical rotational inertia (NCRI) data obtained in torsional oscillator 
experiments^!^ has been unravelled. Specifically, the value of the NCRI and shear modulus 
of helium increases about 2 % as the temperature is lowered down to 0.01 K. Despite 
Beamish et a/.'s findings have been initially rationalized in terms of pinning (unpinning) of 
dislocations induced by the presence of static (mobile) ^He impurities,-*^"— it remains to be 
clarified whether the cited experimental similarities must be regarded simply as coincidental 
or contrarily correspond to manifestations of a same and unique quantum phenomenon 
known as supersolidity.— "— 

Simulation techniques have been demonstrated as invaluable tools for predicting and 
accurately characterizing the energetic and structural properties of quantum solids.—"— 
Nevertheless, computational studies on the elasticity of solid ^He so far have been very 
infrequent. To the best of our knowledgement, there exists only one recent work by Pes- 
soa et al. in which the shear modulus of solid helium has been explicitly calculated from 
first-principles.— This computational scarcity strongly contrasts with research done in other 
fields like classical solid state theory and/or high pressure physics, where estimation of the 
elastic properties of materials (e.g. strain-stress tensor, Griineisen parameters, vibrational 
phonon frequencies, etc.) is a standard.—"— The likely explanation for such a constrast 
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(besides no particular enthusiasm in the matter prior to Beamish et al. experiments) are 
the technical difficulties encountered in modeling of bosonic quantum effects, namely atomic 
delocalization, anharmonicity and particle exchanges. These quantum atomic effects are in- 
deed crucial to comprehend the nature of solid ^He at low temperatures and, as a matter 
of fact, customary harmonic and quasi-harmonic computational approaches'^"— can not be 
used to obtain a reliable picture of it. 

Here, we present a computational study of the elastic properties of perfect (e.g. free 
of defects) solid ^He in the hep structure based on the diffusion Monte Carlo approach. 
This work is intended to improve our understanding of how solid helium reacts to external 
strains/stresses, and further extends the work initiated by Pessoa et alM In particular, we 
provide the zero-temperature dependence of helium elastic constants and related quantities 
(e.g. sound velocities, Griineisen parameters and the Debye temperature) on pressure up to 
~ 110 bar. Our results are compared to experimental data and other calculations whenever 
is possible and, as it will be shown later on, good agreement among them is generally found. 
The computational method that we employ is fully quantum and virtually exact (i.e., in 
principle only affected by statistical uncertainties) so that from a technical point of view 
our study also represents an improvement with respect to previous ffist-principles work- 
based on variational Monte Carlo calculations (i. e., subjected to statistical and importance 
sampling biases). 

The remainder of the article is as follows. In Section [Tll we review the basics of elasticity 
in hep crystals and provide the details of our computational method and strategy. In the 
following section, we present our results along with some discussions. Finally, we summarize 
the main conclusions obtained and comment on prospective work in Section IIVI 

II. THEORY AND COMPUTATIONAL DETAILS 
A. Elastic constants 

For small strains, the zero-temperature energy of a crystal can be expressed as 
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FIG. 1: (a) Representation of the hep unit cell with primitive translational vectors ai, a2 and as, 
and two-atoms basis set (see text), (b) Sketch of the 200-atoms supercell used in the pure shear 
calculations, built by replicating the hep unit cell 5x5x4 times along the primitive translational 
vectors ai, a2 and as, respectively. 

where Vq and Eq are the volume and internal energy of the undistorted solid, {Cjj} the 
elastic constants and {sj} the strain components defined such that Si, S2 and S3 are fractional 
increases in the x, y and z directed axes, and S4, S5 and sq angular increases of the xy, xz 
and yz angles.—"— The symmetry properties of the crystal under consideration define the 
number of elastic constants which are non-zero. For hep crystals, only five elastic contants 
are different from zero, namely Cn, C12, C33, C13 and C44, where C44 is commonly known 
as the shear modulus and abbreviated fi. In order to calculate these five non-zero elastic 
constants is necessary to compute the second derivative of the internal energy of the crystal 
with respect to the strain tensor aij. For this, the hep crystal must be first considered 
in full symmetry, that is, expressed in terms of its unit cell with primitive translational 
vectors ai = a + -^j j , = a + "^j j and = ck (where a and c are the lattice 
parameters in the basal plane and along the z axis respectively, and i, j and k correspond 
to the usual unitary Cartesian vectors), and two-atom basis set ri = ^ai + |a2 + |a3 and 



ra = (0,0,0) (see Fig.^. 
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The relationships between the hep elastic constants {Cij} and applied strain were deter- 
mined long-time ago within the framework of elasticity theory, those being^^i""— 

K = -V = ^33 (Cii + gi2) - 2Cf, 

ydVjy^y^ Cn + Ci2 + 2C33 - 4Ci3 ' ^ ^ 

f dlnc/a \ ^ C33 - Cn - Cu + C13 . . 

\ dV J y^y^ Cn + Cu + 2C33 - 4Ci3 ' ^ ^ 

Co = Cn + C12 + 2C33 - 4Ci3 , (4) 



and 



C44 — C44 . (6) 

Equations ([2]) and ([3]) correspond to homogeneous strains that change the volume and 
shape of the hep unit cell. The dependence of pressure P and c/a ratio on volume can be 
readily obtained from standard equation of state calculations. On the other hand, equa- 
tions @, dS]) and ([6]) represent heterogeneous strains that keep the volume of the hep unit 
cell unaltered. In order to calculate the value of the pure shears Co, Cqq and C44 is necessary 
to compute the variation of the internal energy of the equilibrium structures with respect to 
certain deformations, which can be written as transformed primitive translational vectors. 
In the Co case, these are 




(7) 

where e = (1 + riY^"^, Cq = (^j^) equilibrium condition is satisfied at r] = 0. 

For Cee, we have 
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ck , (8) 



where Cqq = (^^^) equilibrium condition is satisfied at 7 = 1. And finally 

for C44, 



^2 



44 




= ck , (9) 

where C44 = ^ (^^^) equilibrium condition is satisfied at = 0. 

Once the value of the bulk modulus K and quantities din {c/a)/dV, Cq, Cqq and C44 is 
determined, one can calculate the five corresponding Cij 7^ hep elastic constants straight- 
forwardly by solving the non-linear system of equations defined by Eqs. ([2])-([6]). 



B. Diffusion Monte Carlo 



The fundamentals of the diffusion Monte Carlo (DMC) method have been reviewed with 
detail in other works^i^^ so for brevity's sake we recall here only the essential ideas. 

In the DMC approach, the time-dependent Schrodinger equation of a quantum system of 

interacting particles is solved stochastically by simulating the time evolution of the Green's 
function propagator e~^^* in imaginary time r. For r — i- 00, sets of configurations (walkers) 
{Rj = ri,...,rAr} rendering the probability distribution function (\l'o^) are generated, 
where \l/o is the true ground-state wave function of the system and \E' the trial wave function 
used for importance sampling. Within DMC, virtually exact results (i.e., subjected to 
statistical uncertainties only) are obtained for the ground-state total energy and related 
quantities in bosonic quantum systems.—"— It is worth noticing that despite asymptotic 
DMC values do not depend on the choice of the trial wave function, the corresponding 
algorithmic efficiency is tremendously affected by the quality of \E'. 

We are interested in studying the ground-state of perfect hep ^He, which we assume to 
be governed by the Hamiltonian H = —2^^ Sili + Si<j ^e-He('^ij) where is the 
mass of an helium atom, the number of particles and Vne-He the semi-empirical pairwise 
potential due to Aziz et al.^'^ It is worth noting that this two-body potential provides an 
excellent description of the He-He interactions, including weak long-ranged van der Waals 
forces, over all the pressure range considered in this work.— 
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The trial wave function that we use for importance samphng \1/snj simultaneously repro- 
duces crystal ordering and Bose-Einstein symmetry (that is, remains unchanged under the 



permutation of atoms). This model wave function was recently introduced in Ref. j39| and 

it reads 

N N / N \ 

=n/Mn ' (10) 

i<j J=l \t=l / 

where the index in the second productory runs over perfect lattice position vectors (sites). 
In previous works, we have demonstrated that \1/snj provides an excellent description of the 
ground-state properties of bulk hep ^He^ and quantum solid films.— ^^i^ The key ingredient 
for this progress stays in the \1/snj localization factor (second term in Eq. f IlOp ). which is 
constructed in such a way that voids originated by multiple occupancy of a same site are 
energetically penalized. Correlation functions in Eq. f lTOj) were adopted in the McMillan, 
/(r) = exp[— 1/2 (b/r)^] , and Gaussian, g{r) = exp[— 1/2 (ar^)], forms. The value of 
the parameters in factors / and g were optimized variationally at density p = 0.480 
{a = 2.556 A, b = 1.08 a and a = 10.10 cr~^) and kept fixed in the rest of simulations 
performed at different densities. 

The technical parameters in our calculations were set in order to ensure convergence of 
the total energy per particle to less than 0.02 K/atom. For instance, the value of the mean 
population of walkers was held to 400 and the length of the imaginary time-step Ar was 
5 ■ 10~^ . Statistics were accumulated over 10^ DMC steps performed after system 
equilibration, and the approximation used for the short-time Green's function e" -f.^^- was 
accurate up to order (Ar)^.— 

C. Computational strategy 

In order to work out Eq. (j2]) we used the bulk modulus volume dependence reported in 
Ref. {42], where the equation of state of perfect hep '^He was already calculated employing 
the DMC method and considering variational finite-size corrections to the total energy.— 
The equilibrium value of the c/a ratio was found to be constant and equal to 1.63(1) over 
all the pressure range < P < 110 bar. (This outcome is consistent with previous first- 
principles results obtained by other authors.—) Consequently, the left-hand side of Eq. 
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FIG. 2: C44 shear energy results obtained in perfect hep ^He at density p = 0.480 a"^. The 
equihbrium value corresponding to the undistorted hep structure is quoted at e = and the solid 
line represents a third-order polynomial fit to the energies. 

vanishes and the solution to the fifth-order equation system defined by Eqs. ([2])- ([6]) is 



Cii — K + Cge + 7^C*o 

io 

C12 = K — Cee + 7^C*o 

io 

Ci3 = K - -Co 

C33 = K + -Co 
C44 = C44 . 



(11) 



The simulation box used in our pure shear calculations contains 200 He atoms and was 
generated by replicating the hep unit cell 5 times along the ai and a2 directions and 4 
times along the c axis (see Fig. [T]). In proceeding so, hexagonal symmetry in our supercell 
calculations is guaranteed by construction. Periodic boundary conditions were imposed 
along the three directions defining the edges of the non-orthorombic simulation boxes. 

The value of the second derivatives involved in Eqs. (H])-© were computed following 
the next strategy. For each volume and pure shear considered, first we calculated the 
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total energy per particle in a series of supercells generated by incrementally distorting the 
equilibrium geometry according to the transformed translational lattice vectors ([7])- ([9]). Up 
to eight different and equally spaced shear increments (e.g., e, 7 and 0) were considered for 
each volume, taking both positive and negative values. Subsequently, the series of shear- 
dependent total energies so obtained at fixed volume were fitted to a third-order polynomial 
function of the form f{x) = a + bx"^ + cx^. In all the cases, we found that the optimal a, b and 
c values reproduced the series of calculated total energies per particle within their statistical 
errors (see Fig. [2]). Regarding total energy shifts correcting for the finite-size effects, no 
variational energy corrections were considered in the pure shear calculations. The reason 
for this is convenience since finite-size energy corrections are not expected to depend on the 
small shear distortions considered in this work, and consequently they do not contribute to 
the value of the second derivatives involved in the calculation of the elastic constants. We 
shall stress that the volume of the simulation cell in pure shear calculations is kept fixed 
in constrast to bulk modulus and dln{c/a)/dV calculations, in which accurate finite-size 
corrections to the energy are certainly required. 

III. RESULTS AND DISCUSSIONS 
A. Elastic constants 

In Figures [3l H] and [5], we show the pressure dependence of the five elastic constants 
of perfect hep ^He as obtained in our calculations. The error bars 6Cij in our results, 
steming from both the statistical uncertainties of the energies and corresponding third-order 
polynomial fits, typically amount to 6Cij/Cij ~ 2% . We found that the pressure variation 
of all five elastic constants is practically linear within all the studied range. Consequently, 
we performed fits of the form Cij{P) = aij + bijP to our results (see Figures [3l H] and [5]) , and 
obtained a series of Ojj and bij coefficients that we quote in Table I. We note that average 
variance values obtained in the reduced-x^ tests corresponding to our fits were always smaller 
than 2. 

Comparison between our DMC calculations, previous variational Monte Carlo (VMC) 
results and experimental data is provided also in Figures [3H5]. VMC results reported at P ~ 
34 bar-22 have been obtained by Pessoa et al. using a shadow wave function model (SWF).— 
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FIG. 3: Zero-temperature C\\ and C12 elastic constants of perfect hep '^He as a function of pressure. 
Previous variational Monte Carlo (VMC) calculations [20] and experimental data Q] are shown 
for comparison. The vertical dotted line represents the zero-temperature freezing pressure of '^He 
and the straight dashed lines are linear fits to the DMC results (see text). 
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FIG. 4: Zero-temperature C13 and C44 elastic constants of perfect hep '^He as a function of pressure 
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Previous variational Monte Carlo (VMC) calculations 
shown for comparison. The vertical dotted line represents the zero-temperature freezing pressure 
of ^He and the straight dashed lines are linear fits to the DMC results (see text). 
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FIG. 5: Top: Zero-temperature C33 elastic constant of perfect hep ^He as a function of pressure. 



A previous variational Monte Carlo (VMC) calculation 



20l | and experimental data 



are shown 

for comparison. The vertical dotted line represents the zero-temperature freezing pressure of ^He 
and the straight dashed line is a linear fit to the DMC results (see text). Bottom: Dependence of 
the C33 elastic constant on volume. The dashed line represents a power law fit to the DMC results 
from which the value of the corresponding Griineisen parameter is obtained (see text). 
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TABLE I: Value of the parameters obtained in the linear fits to our Cij{P) results (see text), a^j's 
are expressed in units of bar. 

This type of trial wave function correctly accounts for the atomic Bose-Einstein statistics, is 
translationally invariant and so far it has yielded the most accurate variational description 
of solid helium.— Arguably, Pessoa's VMC predictions are in fairly good agreement with 
our DMC results since in general the relation \CY^^ - Cf^^\/Cf^^ < 10% is fulfilled. 
(This inequality is only violated by C12 however in that case measurements appear to follow 
closely our results.) Recalling that evaluation of Ci/s requires from the computation of total 
energy second derivatives, it can be said that the satisfactory DMC- VMC agreement found 
further corroborates the excellent variational quality of the SWF model. 

Regarding the experimental data taken from Refs. 2|j3|, we also find good agreement 
with them (see Figures • The sound- velocity measurements performed by Crepeau and 
Greywall involved high-quality single helium crystals whose basal plane orientations were 
accurately determined using x-rays. Consequently, our modest discrepancies with Crepeau 
and Greywall's data are very likely to be originated by thermal effects given that the temper- 
ature conditions in those experiments were T ~ 1 K. Reassuringly, our C44 results reproduce 
closely recent ^He shear modulus measurements performed by Beamish et al. at just few 
mK (see Figure H]) .-1^ We will comment again on this issue in Section UlI Bl however it can be 
already claimed that the manifested overall good accordance between our Cij calculations 
and 30 < P < 60 bar experiments appears to endorse the rehability of our computational 
approach. 

Another quantity of interest in the study of crystal elasticity is the Griineisen parameter 
7. Essentially, this parameter quantifies how atomic vibrations in a crystal are affected 
by changes in volume. This quantity is customarily defined by 7 = {V jCy) PK , where 
stands for the specific heat and /3 for the coefficient of the thermal expansion. However, 
this definition is not practical for low temperature calculations since in general quantities 
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C„ and P tend to zero in a similar trend near T = 0, and thus leads to an indetermination. 
Alternatively, Klein et alM derived a Griineisen parameter expression that is valid in the 
zero-temperature limit and which depends on the individual vibrational frequency modes. 
Specifically, Klein's expression can be reformulated in terms of the elastic constants as^ 

2d\nV 6' ^ ^ 

Naturally, the volume dependence of each Cij elastic constant can then be fitted to a function 
of the form 



SO that one can readily obtain the value of the corresponding jij parameter. We proceeded 
in this way using the Ctj results obtained in our < P < 110 bar simulations (see Figure [5]) 
and got 7n = 2.34 (5) , 712 = 3.69 (5) , 713 = 4.23 (5) , 733 = 2.70 (5) and 744 = 1-91 (5) , 
where the numerical uncertainties are expressed within parentheses. The averaged Griineisen 
parameter 7 = | X] 7«i corresponding to our results is 2.67 (5), where the summation runs 
over indexes 11, 12, 33, 44 and 66 (766 = 2.70 as obtained from Cqq) because the respective 
elastic constants are the quantities which are directly measured in sound-velocity experi- 
ments.- Our 7 value compares very well with Greywall's experimental result of 2.7, however 
we note that in our calculations the value of the dispersion | ^ \'~fij — 7I is non-zero. It is 
worth comparing the value of the zero-temperature Griineisen parameter of solid helium to 
that of other rare-gas species. We know from Ref. 0] that 7 is 2.5 in Ne, 2.7 in Ar, 2.7 
in Kr and 2.5 in Xe. Consequently, the elastic constants of all five noble gases will vary 
very similarly upon a same volume change. The same conclusion, however, does not apply 
to pressure since the bulk modulus of each element is appreciably different from that of the 
others. 

In order to quantify the importance of quantum effects in our study, we computed the 
contribution of the potential and kinetic energies to the shear modulus (Cf4 and C44, re- 
spectively). For this, we carried out simulations at density p = 0.480 a"^ in which the 
exact value of the second derivative of the potential energy with respect to strain was 
calculated using the pure estimator technique.— The kinetic energy contribution to the shear 
modulus C44 was subsequently obtained by subtracting the quantity C^^ = ^ (^^^ j 
to C44.(We checked that the strain dependence of could also be accurately fitted to a 
third-order polynomial function). In fact, the T = value of C44 in a classical crystal exactly 
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FIG. 6: Pressure dependence of the longitudinal (L) and transvere (T) sound velocities of hep ^He 
along its corresponding c-axis and basal plane, vd represents the averaged Debye velocity (see 
text). Basal sound- velocity data from Refs. Q| (A), (a) and {ij] (v) are shown for comparison. 
The vertical dotted line represents the zero-temperature freezing pressure of ^He. 

amounts to zero since the atoms there remain frozen in their perfect lattice positions (that 
is, = C44). Even in the case of considering quasi-harmonic zero-point motion corrections 
to C44, C44 is not expected to depart significantly from zero. In contrast, we found that in 
perfect hep ^He Cl^^/C^^ amounts to 68 %, or conversely, C44/C44 = 32 %. This last result 
quantifies the quantum nature of solid helium's elasticity and demonstrates the inability of 
classical and quasi-harmonic approaches for reproducing it. 



B. Sound velocities 

Sound velocities in solids, either longitudinal or tranverse, depend on their direction 
of propagation. In crystals with hexagonal symmetry two main propagation modes are 
identified, one along the c-axis (defined by vector as in Section Hi Al) and the other contained 
within the basal plane (defined by vectors ai and sl2 in Section III A|) . The relationships 
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FIG. 7: T = Debye temperature of hep ^He as a function of pressure (solid line). Experimental 
results from Refs. {si] and 52] are shown for comparison. The vertical dotted line represents 
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uncertainty associated to our calculations (e.g. 5Qd/G)d ~ 1%) ■ 



between the elastic constants and sound velocities in hep crystals are^i 

VT2 = (C44/p)^/^ 



48,49 



(14) 



along the c-axis, and 



VL = {Cu/pf 

Vti = {C^Q/pf^ 

VT2 = {CM/pf^ 



(15) 



within the basal plane. 

In Fig. [6l we plot the pressure dependence of the tranverse and longitudinal sound ve- 
locities of hep ^He as obtained from our Cij results reported in Section IIII A[ The error 



16 



bars in our results, not shown in the figure, are 5vl,t/vl,t ~ 1% • It is observed that at 
compressions far beyond freezing all sound velocities increase almost linearly with pressure. 
In contrast, the longitudinal c-axis and basal components appear to follow a certain power- 
law within the low-density interval < P < 25 bars. Certainly, the nature of the sound 
propagation modes in metastable solid ^He, either at positive or negative pressures, is poorly 
understood at present in spite of its fundamental physical interest.— It is our aim to report 
with detail on this topic in elsewhere so we leave discussions on this matter out of this work. 

Experimental longitudinal and tranverse basal sound velocities are shown for comparison 
in Fig. O The agreement between those measurements and our predictions is generally good 
(in fact, as good as claimed in the previous section for the elastic constants). Specifically, 
our predicted sound velocities are systematically a bit larger than those values reported by 
Wanner^, Crepeau^ and Greywall.^ Such a systematic overestimation is consistent with our 
previous suggestion that certain thermal effects, non-reproducible by our simulations, could 
be affecting the experiments. As a fact of matter, the less rigid a material becomes by effect 
of temperature, the more slowly sound waves propagate across of it. 

In order to provide a more meaningful comparison between our zero-temperature re- 
sults and experiments, we computed the T = Debye temperature of ^He The zero- 
temperature Gd of a crystal can be easily extrapolated from lattice heat-capacity mea- 
surements performed at low temperatures, and results for this quantity have already been 
reported for helium over a wide pressure range.— The definition of the T = Debye 
temperature is 

27rn / 3 \ ^ 

0^ = ^ ^D, (16) 



ks \A-nV ^ 

where V is the volume per atom and vo the Debye velocity. This velocity is given by 

where the average velocities vl and vt are defined by 

=3^ = (^) ' (18) 

and the (■ ■ ■ ) brackets denote angular averages of the longitudinal and tranverse velocities. 
In our case, we have approximated the angular averages in Eq. [TSlby 

1 / 1 \ 1 1 , 

(19) 
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and 

4 "3 (4;;) +3 (4;;) +3 (4:) ' 
where index b stands for basal and index c for c-axis.— 

In Fig. [71 we plot ou r resu lts for the zero-temperature Qd of hep ^He and experimental 
data taken from Refs. 
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52| . In fact, excellent agreement between Gardner and Ahlers 



measurements and our calculations is observed. This last result appears to further ratify our 
previous suggestion that, once thermal effects are corrected for, our T = elastic constants 
and deriving quantities predictions closely reproduce experiments. Finally, we note that the 
pressure variation of 0d is very similar to that observed in the longitudinal sound- velocity 
components of helium, namely almost linear at high compressions and of power-law type at 
low densities. 



IV. SUMMARY AND PERSPECTIVES 



We have developed a fully quantum computational strategy to accurately calculate the 
zero-temperature elastic constants of perfect hep ^He under pressure. Our diffusion Monte 
Carlo results are shown to be consistent with low-T sound-velocity measurements and previ- 
ous variational first-principles calculations. It is found that all five non-zero elastic constants 
of helium vary linearly with pressure within the range < P < 110 bar, and we have provided 
an accurate parametrization of each of them. The Griineisen parameters, sound velocities 
and T = Debye temperature of solid helium have been also determined and compared to 
experimental data. The computational method introduced in this work is totally general so 
that it can be used for the study of any other hep quantum solid appart from helium (e.g. 
H2), and/or be conveniently altered in order to investigate other crystal structures (e.g. face 
centered and body centered cubic). 

It is our intention to analyze the elastic behavior of "^He at negative pressures using the 
computational technique described here. In doing this, we expect to be able to determine its 
spinodal density limit (that is, the density at which the elastic constants vanish) rigorously, 
and also characterize the pressure dependence of the tranverse and basal sound velocities 
near it. Also we are interested in applying our formalism to the study of the ground-state 
of defective hep "^He (for instance, by introducing vacancies), where the supersolid state of 
matter clearly manifests. In doing this, we expect to gather quantitative knowledge on the 

18 



relationship (if any) between elasticity and supersolidity and so to help to understand the 
origins of recent shear modulus observations. Work in these directions is already in progress. 
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